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Abstract. We present a generalized Keller-Segel model where an arbitrary number of 
chemical compounds react, some of which are produced by a species, and one of which is a 
chemoattractant for the species. To investigate the stability of homogeneous stationary 
states of this generalized model, we consider the eigenvalues of a linearized system. 

■ We are able to reduce this infinite dimensional eigenproblem to a parametrized finite 
. dimensional eigenproblem. By matrix theoretic tools, we then provide easily verifiable 

■ sufficient conditions for destabilizing the homogeneous stationary states. In particular, 
^ . one of the sufficient conditions is that the chemotactic feedback is sufficiently strong. 

Although this mechanism was already known to exist in the original Keller-Segel model, 
here we show that it is more generally applicable by significantly enlarging the class of 
models exhibiting this instability phenomenon which may lead to pattern formation. 



< 



1. Introduction 



■ Pattern formation is a fascinating area of biology with many unanswered questions. 

An early example is furnished by morphogenesis, which focuses on the following question: 
How can an initially spherically symmetric system, such as an embryo in early stages of 
CN I development, lead to a creature as spherically asymmetric as a person, or an animal with a 

spotted coat? One of the first attempts to theoretically explain morphogenesis (by means 
^ . of a mathematical model) is by Alan Turing. In [18] , he proposed that an instability 

^ i mechanism in specific systems of chemical reaction-diffusion systems is the underlying 

cause of morphogenesis by showing that under certain conditions on both the reaction 
and diffusion terms, some such systems which are stable without diffusion are destabilized 
with diffusion [12]. Over time, pattern formation in various biological systems has been 
attributed to such Turing instabilities in systems of reaction-diffusion equations. 

The main goal of this paper is to present an alternate avenue potentially leading to 
^ ■ pattern formation via chemotaxis. The key idea of this mechanism is already present in 

the well-known Keller-Segel model [10] (reviewed below), as elucidated by Schaaf [16] . 
It relies on the destabilization via a chemotactic term of an otherwise stable uniform 
steady state. This mechanism of destabilization based on the chemotactic sensitivity is 
also analyzed for the case of the two equation system in one spatial dimension in [TH] . In 
this paper, we will identify a larger class of systems in higher spatial dimensions where a 
similar chemotactic destabilization mechanism may lead to pattern formation. 
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Chemotaxis occurs when the movement of a species is influenced by chemicals in the 
environment [7]. It naturally arises in a wide variety of interesting biological settings, in- 
cluding cancer metastasis, angiogenesis, immune system function and egg fertihzation [6] . 
Considering reaction diffusion systems with chemotactic terms, we show that there are 
quite general classes of chemical reaction networks (CRNs) which can be destabilized by 
means of an inherent feedback mechanism. The feedback arises as follows. A population 
of cells produces a chemical, which in turn participates in a CRN. The CRN leads to one 
or more chemical products. One of these products serves as the chemotactic signal for 
the cells, thereby completing the feedback loop. Under certain conditions this feedback 
loop has unstable homogeneous steady states which, in the absence of the feedback, would 
be stable. The destabilization of homogeneous stationary states, caused by the feedback 
mechanism built into the system, may give rise to pattern formation. As mentioned ear- 
lier, such a mechanism is already present in a simple CRN appearing in the well-known 
Keller-Segel model |10]. In this model, the signal emitted by the cells is the same as the 
chemotactic signal, and it is subject only to decay. Since this CRN consists of a single 
decay reaction, it is perhaps one of the simplest possible CRNs. In this paper, we find a 
larger class of CRNs with destabilizable homogeneous steady states. 

To further explain the difference with the traditional destabilization mechanism, let us 
first review Turing's diffusion driven instability for a two-species system (generalizations to 
multi-species systems are also well-known), as explained in [13], using a reaction-diffusion 
system of the form 

(la) Ut = diAu + f{u,v), X E Q, t>0 

(lb) Vt = d2Av + g{u, f ), x G fi, t > 

Oil f)l) 

(ic) ^ = :^ = o, xedn, t>{). 

on on 

where di and d2 are positive diffusion coefficients. Here Ut denotes du/dt and n (in 
d/dn = n ■ V) denotes the outward unit normal on the smooth domain boundary dQ. 
Assume the existence of a homogeneous steady state {u*,v*) of ([1]). When the diffusion 
terms are ignored, we obtain the ordinary differential equation system corresponding to 
just the reaction terms: 

(2a) — = f(u,v), t>0 

(2b) ^ = g(u,v), t>0. 

We note that {u*, v*) is also a steady state of ([2]) and we assume that it is linearly stable, 
i.e., we assume that the Jacobian matrix 



J 



gu{u*,v*) g^{u\v*) 



has eigenvalues with negative real part. (Again, derivatives are indicated by subscripts.) 

The question Turing asked is whether («*,!>*), considered as a steady state of ([I]), can 
be linearly unstable, even if it is a stable steady state of ([2]). In other words, is it possible 
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that the eigenvalue problem arising from linearizing ([T]) at {u*,v*), namely 

d,AU + fuiu*,v*)U + Mu*,v*)V = XU 
d2AV + gu{u*,v*)U + g,{u\v*)V = XU 

d£_d\^_ ^ 

dn On ' 

has an eigenvalue A with positive real part? It turns out that a first necessary condition 
for this to happen is that di ^ d2, implying that both species must move with different 
diffusion constants. By an appropriate scaling we may assume, without loss of generality, 
that di = 1, and then a second necessary condition is that J has one of the sign patterns 



(3) 
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with additional restrictions on d2 and fi. The relevance of this sign pattern is that it 
imposes certain restrictions on how the species interact. For example, the first matrix 
implies that the first species activates itself and the second species, and that the sec- 
ond species inhibits itself, as well as the first species. A generalization of Turing-type 
instabilites to a system modeling an arbitrary number of chemicals can be found in [15j . 

The model that we propose for pattern formation on the other hand, is a fully nonlinear 
react ion- diffusion system which also contains a nonlinear chemotaxis term. Destabiliza- 
tion of homogeneous steady states in our model does not require any restrictions on the 
diffusion coefficients, although we do impose some conditions on the sign pattern of the 
Jacobian corresponding to some of the reaction terms. These may however be less crucial 
than the role played by the chemotactic term. 

The Keller-Segel model for chemotaxis is one of the most widely studied models. It was 
developed focusing on the cellular slime mold Dictyostelium discoideum. As explained in 
[7] and [To], the derivation of the Keller-Segel model is based on the aggregation stage 
of D. discoideum''^ life cycle. The unicellular slime mold grows by cell division until it 
depletes its food source, and begins to enter a starvation mode. After some time, one 
cell will emit a signal of cyclic Adenosine Monophosphate (cAMP). The other cells are 
chemotactically attracted to this signal, and in turn begin to emit cAMP themselves. 
The cells also produce an enzyme which degrades the cAMP by first binding to it and 
forming a complex which in turn breaks up into the enzyme plus what we call a "degraded 
product". To model this, the following assumptions are made in [7] and |10j : 

• Let u{x,t) denote the cell density of the slime mold, v{x,t) the density of the 
chemoattractant cAMP, rj{x,t) the density of the enzyme, c{x,t) the density of 
the complex formed by the binding of cAMP and the enzyme, and d{x, t) denote 
the density of the degraded product. 

• The chemoattractant is produced by the amoeba at a rate f{v) per amoeba and 
the enzyme is produced at a rate g{v,r]) per amoeba, which, according to [10] is 
allowed to vary with both the chemoattractant and enzyme concentrations. 

• The chemoattracant, enzyme and complex react as shown below. 

V + T] i )• C > T] + d. 

We assume that the reaction rates proceed according to the law of mass action. 
Accordingly, let the positive rate constants be denoted by ri (for the reaction 
V + rj ^ c), r^i (for the reaction c — )■ f + 77) and r2 (for the reaction c ^ rj + d). 
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• The cAMP, the enzyme and the complex all diffuse according to standard Pick's 
law. 

• The cell concentration changes by random diffusion, as well by chemotaxis, due 
to the cells "climbing the gradient" of the chemoattractant cAMP. 

The last assumption is again modeled using Pick's law of diffusion, but with the amoebic 
flux J*^") consisting of a purely diffusive part and a chemotactic contribution, namely 

for some chemotactic sensitivity function k2{u,v) and diffusion constant ki. These as- 
sumptions lead to the full Keller-Segel model [7] on a bounded domain Q for time t > 0: 



(4a) Ut = V ■ {kiVu - k2{u,v)Vv), x e Q, t>0 

(4b) Vt = k^Av — rivr] + r^ic + uf{v), x E Q, t>0 

(4c) rjt = kr^Ar] — riVT] + {r_i + r2)c + ug{v,ri), x E Q, t>0 

(4d) Ct = kcAc + riVT] — (r_i + r2)c, x E Q, t > 

, , , du dv dn dc 

4e = =_^= =0, xedn, t>0, 

on on on on 

Here, k^, krf and kc denote the positive diffusion constants for cAMP, enzyme and complex 
respectively. 

In [IQl Keller and Segel use a steady state assumption to reduce the above set of four 
equations to a two-equation system which models only the density of the chemoattractant 
and the cell density. Purther simplifications result in the so-called [6l [7j "minimal system" , 



(5a) Ut = V ■ iyu - xuVv), x e Q, t>0 

(5b) Vt = k^Av — ■yv + au, x E fl, t > 
Ou dv 

(5c) _ = _ = 0, xedn, t>0 
on on 

(5d) u{0,x) = uq{x), v{0,x) = vo{x), x E Q. 

This system has been studied extensively, and there are many results on local and global 
existence, positivity, and blow-up of solutions as well as steady states [IH El [71 [T6] . 
Solutions that exhibit blow-up, defined as solutions for which the L°°-norm of m or v 
becomes unbounded in either finite or infinite time, have attracted considerable interest 
and are reviewed extensively in [7]. (The blow-up phenomenon, restricted to a simplified 
Keller-Segel model, is also featured in [HI Ch. 5], which nicely illustrates the difficulties 
encountered when dealing with this complex issue in an elementary setup.) One significant 
result is that in two dimensions, if ax J^Uo{x)dx < inky then the solution to (jS]) exists 
globally in time and its L°°-norm is uniformly bounded for all time [3 Table 5]. On the 
other hand if Ank^ < ax Jq uo{x)dx < Snk^, then there exist initial data (mq, vq) for which 
the solution of (I5l) blows up at the boundary of Q in finite or infinite time [7| Table 5]. 
Although the question of blow-up of solutions for our generalization of the Keller-Segel 
model is interesting, we do not address that issue in this paper. Instead, here we focus 
exclusively on destabilization of homogeneous steady states. 
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Our work fits in the context of tlie more recent attempts to investigate multi-species, 
multi-cliemical systems witli cfiemotactic terms [21 IH El [IZl EHl El IH] ■ Just as in tlie case of 
tlie Keller-Segel model, the blow-up phenomenon in these generalized models has received 
considerable attention [201 [3l HI IH] • Of particular interest to the issue we address here, is the 
recent work by Fasano et al [5], Horstmann [8], and Liu et al [11]. All papers consider the 
existence of non-homogeneous or periodic steady states, a signature of pattern formation 
for multi species models. In [5] , a two-species two-chemical system is proposed where each 
species produces one of the two chemicals, which in turn serve as attractive chemotactic 
signals for the species that produce them. In addition however, these chemicals are also 
chemotactic repellants for the species that do not produce them. The model is analyzed for 
scalar spatial domains, and focuses on the existence of nonuniform and periodic steady 
state patterns. Motivated in part by [S], Horstmann in [S] proposed an n species, m 
chemicals model where all chemicals are allowed to serve as chemotactic signals for the 
species, and the species themselves can chemotactically attract or repel each other. In 
addition, the chemicals may be produced and/or consumed by the species, and in principle, 
could react with each other in fairly general ways. Several natural questions that have 
attracted the attention of several researchers in the context of the Keller-Segel model, 
are addressed in [8j for the multi-species, multi-chemical models. They include blow-up, 
global existence of solutions, existence and non-existence of nonuniform steady states, 
destabilization of uniform steady states as a way to achieve pattern formation, and the 
existence of Lyapunov functionals. Although [S] begins with this very general model, 
the focus of the paper quickly turns to specific cases, motivated by the literature, where 
only a limited number of non-reacting chemicals (often just one or two) appear in the 
model. One notable exception is in Section 3.1 of [8], where a one-species, two chemical 
system is proposed, where the two chemicals react in a nontrivial way. The model is 
inspired by work of [2] where an E. coli population is studied in the presence of two 
chemotactic attractants, namely glucose and oxygen. In [llj a system with two chemical 
agents is also considered, but in this case one is a chemoattractant while the other acts 
as a chemorepellent. In contrast, the emphasis of the work presented here, is on possible 
pattern formation in models with several chemicals, which react according to fairly general 
reaction networks, but only a single species that chemotactically responds to one of these 
chemicals. 

In Section [2] of this article, we explain our generalized model and the system of partial 
differential equations which arises from it. In Section |3l we discuss conditions for the 
existence of steady state solutions to the system in a special case. In Section H] we state 
our main results, which give sufficient conditions under which steady states are unstable. 
Section El gives examples of theoretical biological systems which meet the conditions of 
our theorems and also provides a motivating example for generalizing our results. 



2. The generalized model 

In this section we present a generalization of the Keller-Segel model (jSj) where the 
modeled species interacts with several chemical compounds. 

As before, consider a population of some species occupying Q, an open bounded con- 
nected ra- dimensional domain. Let the population density at a point xinfl and time t be 
u{x,t). The population climbs the gradient of a chemical of concentration vn- To model 
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the more general and the more hkely biological scenario, where this chemical is in reac- 
tion with several other compounds in the environment, we introduce the concentrations 
of — 1 other chemicals, denoted by fi, f2, • • • , vn-i (see Figure [I]). These chemicals can 
interact with each other as well as with vn creating a CRN. Writing t? = [t>i, . . . , f at]*, we 
model the rate of change of v due to the CRN by g{v) for some function g : — )■ M^. 
In applications, this function g is often determined by mass action kinetics. Note that in 
general g includes decay of the chemicals. 

Some of these compounds can also be produced by the species. To model this scenario, 
consider a general subset 5* C {1,2, ...,A^} and assume that the species produce the 
chemical vi at a rate a; > for all / G S. The remaining compounds are not generated 
by the species, so we set a; = for / ^ S. Let a e be the vector whose ith 
component is the introduced above. Thus, the rate of change of v would be governed 
by Vt = du + g{v), in the absence of any diffusive effects. 




Figure 1. A schematic of the species {u) and the chemicals (fj's) interacting. 

However, we also want to account for spatial diffusion. In general, both the species 
population and the chemical concentrations diffuse. Assume that the diffusion coefficient 
for u is given by the constant D > 0, and the diffusion coefficients for each Vi are given 
by the constants Di > 0. Let D be the diagonal matrix whose ith diagonal entry is Di. 
Assuming that the chemotactic sensitivity function is X'^? for some constant x > 0, we 
arrive at our generalized model: 

(6a) Ut = V ■ {DVu - xuVvn) x eVL, t > 0, 

(6b) Vt = DAv + du + g{v) x eVt, t > 0, 

(6c) ^(x, 0) = 'Uo(x), v{x,Q)=Vo{x) X G fi, 

(6d) dui^^dv^^^ ^^^^^ ^^^^ 1<^<N. 

on on 



Here we have assumed zero Neumann boundary conditions. The initial conditions uq and 
f are assumed to be nonnegative functions on Vt. Notice that upon integrating ( l6a|) over 
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the domain and applying the divergence theorem and boundary conditions, we have 



/ utdx = I V ■ {DVu — xuVvn) dx 
Jn Jn 



{DVu — xuVvn) ■ n dS 



an 
= 0. 

Hence if u is sufficiently smooth we may interchange the integral with the differentiation 
with respect to time to conclude that the total population of the amoeba inside the domain 
is constant in time. 

As an immediate example of an application of our generalized model, consider the 
unreduced Keller-Segel model (jl]) with the chemotactic sensitivity function k2{u,v) = xu 
and the rates f{v) and giv^rf) assumed to be constant. We change the notation to let Vi 
denote the density of the emitted enzyme, the density of the chemoattractant cAMP 
and V2 the density of complex formed from the cAMP and the enzyme. Then setting 



-riViV^i + + r2)f 2 
riviv^ - (r_i + r2)v2 



and 



a 



Oil 



as 



for constant rates > 0, we recover (jlj) in the framework of our model. 

In the following sections, we will examine the constant stationary states of the model ([6]). 
We will be especially concerned with conditions for the destabilization of such constant 
states. 



3. Existence of homogeneous steady states 



A steady state 
(7a) 
(7b) 

(7c) 



^u*, tT*) of our generalized model 
DAu^ — V ■ (xw*V['U*]Ar) : 
DAv^ + du^: + g{v^) ■- 

dn On 



satisfies 

X eVL, 
X ^VL, 



0, 
0, 



— ^ = 0, xedVt. 



If the steady state is homogeneous, i.e., if it consists of spatially constant functions, 
then (ITcI) obviously holds, and ( I7al) -(f7bl) reduces to 

(8) du^ + g{v^) = 0. 

Recalling that the total population of u is conserved, we note that the initial conditions 
will predetermine the value of u^,. A basic question is whether such nonnegative constant 
steady states (u*, -y*) exist. In this section, we will answer this question in the affirmative, 
assuming a certain linearity condition on g. 

Before we do so, we need to introduce some terminology from matrix theory that we 
will use throughout. For any square matrix A, let p{A) denote the spectral radius of A. 
We write A> for a real matrix (or vector) A if each entry of A is nonnegative. Similarly 
A > signifies that every entry of A is positive. 
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Definition 3.1. A Metzler matrix A G M^^^ is a matrix liaving nonnegative off-diagonal 
entries, i.e., 

Aij > 0, « ^ J. 

We will also need M-matrices, which can be defined in many equivalent ways. Follow- 
ing [T], we introduce the definition below. 

Definition 3.2. A matrix A E x is an M-matrix if it satisfies 

<0, i^j 

(2) Au > 

(3) There exists a nonnegative matrix B > and a number s > p{B) such that 
A = sl -B. 

Clearly, if A is an M-matrix, then —A is a Metzler matrix. An M-matrix A is invertible 
if and only if there exists a S > and s > p{B) such that A = si — B. This is because 
by the Perron- Frobenius Theorem [T], the spectral radius of a matrix i? > is also an 
eigenvalue of the matrix. We will also need the following well-known result [B, PP- 137]: 

Lemma 3.3. Suppose A G R^^^ is nonsingular. Then A is an M-matrix if and only 
ifA-^ > 0. 

We are now ready to address the above mentioned question of existence of homogeneous 
steady states, under the assumption that each component of ^ is a linear function. 

Proposition 3.4. Suppose g{v) = —Av for some A in R^^^ and a > 0. If A is a non- 
singular M-matrix, then for every positive constant u^, there exists a unique nonnegative 
homogeneous stationary solution (m^,, v^,) of ([6]). Conversely, A is a nonsingular M -matrix 
if has a nonnegative homogeneous stationary solution {u^,v^) for every u^, > and 
every a > 0. 

Proof. First, suppose that A is a nonsingular M-matrix. Let be any positive constant. 
Then setting 

(9) = u^A~^d, 

we observe that satisfies ([8]). Hence {u^.,v^) satisfies ([7j). Furthermore, iT* > by 
virtue of Lemma [3.3[ so is a nonnegative homogeneous stationary solution. For 

the given m^,, the component is uniquely determined. Indeed, if (m*,^^) were another 
homogeneous stationary solution, then by ([H]), Ai^^^ = u^a. But by the invertibility of A 
and dnD, this implies v'^ = t7*. This proves the first assertion of the proposition. 

Conversely, suppose {u^,vj^'^) is a nonnegative homogeneous stationary solution ob- 
tained using M^, = 1 and a = ei, the unit vector in the ith coordinate direction. Then, 
by ([8]), Av^^^ = Cj. Letting C denote the matrix whose zth column is v^^\ this implies 
that AC = I, the identity. Hence A is invertible, and moreover, since wi*^ > 0, we have 
A~^ = C > 0. Hence by Lemma [3.31 A is a nonsingular M-matrix. □ 

Having given a sufficient condition under which plenty of nonnegative homogeneous 
steady states exist, we now proceed to study the stability of such states, assuming they 
exist. Note that we no longer assume that g is linear in the ensuing analysis. 
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4. Instability of homogeneous steady states 



In order to analyze the stability of a given homogeneous steady state (u*, -u*), we begin 



by linearizing the spatial partial differential operator in 
eigenvalue of the linearized operator, i.e., A satisfies 



about 



Let A be an 



(10a) 
(10b) 

(10c) 



DAU - xu*^Vn 
DAV + aU + JV 
dU _ dVk 
On On 



XU, 
XV, 

0, 



X e^l, t > 
X eQ, t>0 

xedn, A; = 1,2,..., AT. 



for some nontrivial pair of functions {U{x),V{x)) on Q. Here the N x N (Jacobian) 
matrix J denotes the Frechet derivative of g with respect to v 

dg_ 

dv 



:ir 



J 



i.e., Jki = dgk/dvi evaluated aX v = v^. The purpose of this section is to study the 
eigenvalue problem flTU]) and thereby find conditions under which the homogeneous state 
(m^,, v^) is unstable. 

We will use the weak formulation of (fTOj) , where u and the components Vi are all sought 
in the (complex) Sobolev space H^{Q) consisting of square integrable functions whose 
first order distributional derivatives are also square integrable. By definition of the weak 
eigenproblem, the number A is an eigenvalue of the linearized operator if there exists a 
nontrivial (f/, Vi, . . . Vn) in H^{Q)'^'^^ satisfying 



:i2a) 



- / DVU • Vyp + / XM*VV^ ■V^ = X U^, 
'n Jn Jn 

N 



(12b) - [ DkWk -Vip + au ! U^+y^Jki f ViTp = X [ V^Tp, 

Jn Jn Jn Jn 

for k = 1,2, . . . , N, and for all ip in H^{Q). The integrals are with respect to the standard 
Lebesgue measure (omitted). The system f|T2l) can be obtained from the classical formu- 
lation ffTOj) by multiplying the equations of ffTOj) by a test function Ip and integrating by 
parts, provided the boundary of the domain dfl and the eigenfunctions u, v are smooth 
enough. So as to admit irregular domains with possibly nonsmooth eigenfunctions, we 
adopt the weak formulation ( TT2l) as the definition of our eigenproblem. 

4.1. Reduction to finite dimensions. The first step in our study of ffTOj) proceeds by 
generalizing a method of Schaaf [16] . This allows us to relate the eigenvalues of the infinite 
dimensional eigenproblem f lT2|) to a finite dimensional eigenproblem. The latter is easier 
to analyze. The finite dimensional eigenproblem involves an (A^ + 1) x (A^ + 1) matrix, 
written in block form as 



M(/x) 



DjJ, K 

a J + jiD 

where = (0, 0, ■ ■ ■ , — xm*/^) and /i is a real parameter. The parameter will always be 
set to one of the numbers 



= /io > /il > /U2 > /U3 > 
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in the spectrum of the Laplace operator with Neumann boundary conditions, i.e., Hi 
satisfies 



(13) 



/ VWj ■ V<y9 = /ij / Uiip, V(y9 G H^{Q), i = 0, 1, 2, . . . 
Jn Jq 



where Ui is the corresponding eigenfunction Ui in H^{Q), normahzed to have unit L^(f2)- 
norm. With these notations, the reduction to finite dimension is achieved by the following 
result. 

Theorem 4.1. The number X solves the eigenproblem ( |T2i) if and only if it is an eigenvalue 
of M{fii) for some integer i > 0. 

Proof. Suppose A satisfies (IT^ with a nontrivial {u,Vi, . . .v^) in H^{Q)'^~^^ . Then, since 
{ui : i > 0} form a complete orthonormal set in L'^{Q), it is possible to find an i such 
that not all of the N + 1 numbers 

Xq = McUj, and Xj = / VjUi, for j = 1, 2, . . . , A^, 
Jn Jn 

are zero, i.e., the vector x = [xq,Xi, ■ ■ ■ ,xj\fY is nontrivial. Furthermore, choosing (f = Ui 
(with the above found i) in (fT2!) . using ( !T3|) . and rewriting the resulting equations in terms 
of the numbers x,-, we find that 



'Dfii 

"1 Jll+f^iDi Ji2 

a2 J21 J22 + i^iD2 



J. 



Nl 



J. 



N2 



JlN 
J2N 

Jnn + l^iD 



N 



X = Xx. 



The matrix above is the same as M{fii). Thus we have shown that A is an eigenvalue of 
M{fj,i) with a; 7^ as its corresponding eigenvector. 
To prove the converse, suppose A satisfies 

M{fii)x = Xx 

for some x ^ and some i. Then, setting 

U = XoUi and Vk = x^uJi 



for A; = 1,2,..., A^, 



we find that the equations of (IT^ hold for any if G H^{Q). Thus, A is an eigenvalue of the 
linearized operator with the (nontrivial) pair {u, v) as its corresponding eigenfunction. □ 

Remark 4.2. Recalling that the mean of u in (l6al) does not vary with t, we may restrict 
ourselves while linearizing to perturbations U with j^U = 0. Suppose we also restrict 
the perturbations in to be spatially nonhomogeneous by imposing j^V^ = 0. Then, 
following the proof of Theorem 14.11 it is easy to see that A is an eigenvalue of M(/ij) for 
some 2 > 1 if and only if it solves the eigenproblem (IT^ with the additional conditions 
f^U = 4 ^fc = for all A; = 1, . . . , N. 
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4.2. A sufficient condition for destabilization. We now focus on conditions under 
which the linearized operator has at least one eigenvalue A in the right half of the complex 
plane. In such a case, the corresponding stationary state is called linearly unstable. In view 
of Theorem 14.11 we only need to study the eigenvalues of the parametrized matrix M(yu). 
Let us begin by reviewing a few useful results from matrix theory that we shall need. The 
following definitions and the next two theorems are well-known - see e.g. [H Theorems 1.4 
and 2.35 in Chapter 2]. 

Definition 4.3. The directed graph G{A) associated to an x matrix A, consists of 
vertices Vi,V2, ...,Viy where an edge leads from Vi to Vj if and only if ^ 0. 

Definition 4.4. A directed graph G is strongly connected if for any ordered pair {Vi, Vj) 
of vertices of G (with i j), there exists a sequence of edges (a path) which leads from 
V to Vj. 

Definition 4.5. A matrix A is called irreducible if G{A) is strongly connected. 

Theorem 4.6. (Perron-Frobenius) Let A > be a square matrix. 

(1) If A > 0, then p{A) is a simple eigenvalue of A, greater than the magnitude of any 
other eigenvalue. 

(2) If A is irreducible, then p{A) is a simple eigenvalue, any eigenvalue of A of the 
same modulus is also simple, A has a positive eigenvector x corresponding to p{A), 
and any nonnegative eigenvector of A is a multiple of x. 

Theorem 4.7. (Spectral radius bounds) Let A > be an irreducible N x N matrix. 
Letting Sj denote the sum of the elements of the ith row of A, define 

S(A) = max Si and s(A) = min Sj. 

l<j<Af l<j<Af 

Then the spectral radius p{A) satisfies 

s{A) < p{A) < S{A). 

With the help of these well-known results, we can now give a sufficient condition for 
the destabilization. 

Theorem 4.8. Assume that the system ([6]) has a positive homogeneous steady state solu- 
tion (m*, and let J be as in flTT]) . Suppose J and a > satisfy the following conditions: 

(1) There is an with 1 < i^ < N such that > 0, 

(2) J is irreducible, and 

(3) J IS Metzler. 

Then, if either xu^. or is sufficiently large, then (m*,'U*) is linearly unstable. 

Remark 4.9. The biological interpretation of the first condition is that the organism must 
produce at least one of the chemicals in the reaction network. The second condition is 
a technical condition which we will relax in the next theorem. The third condition is 
the most restrictive of the conditions, but we give an example of a large class of CRNs 
for which it holds in Corollary 15. 2[ We also give examples of CRNs which do not meet 
condition three but for which a similar result still holds (see Section [S]). 
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Remark 4.10. For any fixed i^, we may rescale by u'^ = ai^u^, to rewrite fl7a|) as 

DA<-V-(CV[t;*]7v) = 

where C = ai^x^*- Then the condition on x^* or in Theorem 14.81 can be interpreted 
in terms of the single parameter C . 

Proof. We use Theorem 14.11 Accordingly, it is enough to prove that M{^i) has a positive 
eigenvalue A for some /Zj. We will now show that for any /x < 0, the matrix M(/i) has a 
positive eigenvalue under the given assumptions. 

To this end, fix some < 0, and set K = —fiu^^x > 0- Let 

B = M(/i) + r J, 

where r > is chosen large enough so that B has positive diagonal entries. Note that 
while both r and K are dependent on /i, r does not depend directly on the value of 
K (which may be changed by varying parameters other than fi). Note also that r is 
independent of a. Then since J is Metzler, i? is a nonnegative matrix. 

In fact, B is an irreducible nonnegative matrix. To see this, consider the directed graph 
of J made of vertices Vi, . . . , Vjy. The graph G{B) is obtained after augmenting G{J) by 
a vertex corresponding to the first row and column, say Vq, and adding vertex loops as 
needed to account for the positive diagonal of B. Since J is irreducible, G{J) is strongly 
connected. Moreover, since a 7^ 0, we conclude that there is a path from Vj to Vq for any 
j. Similarly, since K > 0, there is a path from Vq to Vj for any j. Hence, G{B) is strongly 
connected and B is irreducible. 

We now claim that 

(14) lim p{B) = 00 and lim p{B) = 00. 

The theorem follows from the claim. Indeed, since 5 is a nonnegative irreducible matrix, 
p{B) is an eigenvalue of B by Theorem 14.61 Hence, by ([HI), for large enough K, we can 
ensure that p{B) > r. But then, p{B) — r is a positive eigenvalue of M(/i) = B — rl. 
By the same reasoning, for large enough Oj^, the matrix M(/i) has a positive eigenvalue. 
Thus, in either case, the stationary state (m*,^*) is unstable. 

In the remainder of this proof, we establish the claim f|T^ . To this end, let us denote 
by -B* the transpose of the matrix B and note that for any integer m > 1, 

p{BT = piiBT)- 

Therefore, to prove that the first limit in f|T4l) is infinite, it suffices to prove that 

(15) lim p((S*)™) = 00 

for some integer m > 1, take the mth root, and use p{B) = p{B^). 

To prove (1151) . we start by observing that the (i, j)th entry of (5*)*" is 

(16) [{BTh = E E ■ ■ ■ E [^*]^- [^']-.- ■ ■ ■ [^*]--.^- 

n i2 im-l 

where the sums run over the ranges of the matrix indices ii, . . -im-i- We choose these 
ranges to be 0, 1, ... (instead of the customary 1, 2, . . . , + 1) so as to match the 
indices of the vertices Vq, . . . , V)v of G{B). Since G{B^) is strongly connected, there is a 
path from Vi to Vj for any i and j of some length (number of connecting edges) I. Since 
5* has positive diagonal entries, each vertex of G{B^) has a loop, and consequently, if 
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there is a path of length I, then there is a path of any length longer than / as well. Hence, 
there is a number m such that every two vertices in G{B^) are connected by a path of 
length m. With this m, consider the {i, 0)th entry of (S*)™, as given by ( 1T6|) . Since there 
is a path connecting Vi to Vq, say Vi —> Vi^ —> Vi^ - ■ ■ Vi^^_^ Vq with im-i = N, we have 
from (fT6!) that 

[(5*)'"]^o > ■ ■ ■ Bo,N = aiK)K, 

where Ci{K) = Bi^^iBi^ i^ ■ ■ ■ Bi^_^^i^_^ is a positive and nondecreasing function of K (in 
fact a polynomial in K). Let C{K) denote the minimum of such Ci{K) for all < z < A^. 
Then, the minimal row sum s{B) satisfies 

siiBT) > ^^njBTU > C{K)K, 

so by Theorem 14. 7[ 

> > C{K)K. 

Since C{K) is a nondecreasing positive function of letting — oo, we prove ( TT^ . 

This proves that the first limit in f ll4p is infinity. The proof that the second limit is 
also infinity follows in the same manner by finding, for each < i < A^, a path of length 
m from to l^*. □ 

4.3. Relaxing the irreducibility condition. As in the previous proof, let Vi, V2, . . . Vn 

be the vertices of the G{J). One way to check the second condition (irreducibility) of 
Theorem 14.81 is to verify that there is a path from any Vi to Vj for all i ^ j. In this 
subsection, we will give another condition involving only one path that is often easier to 
check. 

Let us first recall some further terminology. A strongly connected component of a graph 
is a maximal strongly connected subgraph. Note that a strongly connected graph consists 
of exactly one strongly connected component. We consider the following equivalence 
relation on vertices of a graph: Two vertices Vi and Vj in a graph are equivalent if there 
is a directed path from Vi to Vj and a path from Vj to V.. Then, the equivalence classes 
created by this equivalence relation consist exactly of the vertices of the strongly connected 
components of the graph. Keeping these in mind, recall the following definition of pQ . 

Definition 4.11. The classes of an N x N nonnegative matrix A are the disjoint subsets 
{ii,i2, • • • , h} ^ {1; 2, ... , A^} corresponding to the vertices of the equivalence classes 
of G{A). We also identify a class with its strongly connected component. 

If A is reducible, then by reordering the vertices, a permutation matrix P can be found 
so that 

(17) T = PAP^ 

is a lower block triangular matrix with the blocks consisting of the classes of A (see pQ). 
Since the classes of A are strongly connected, the diagonal blocks of T are irreducible. 
Using these facts we can drop the irreducibility assumption in Theorem 14.81 and replace 
it with a weaker assumption as follows. 

Theorem 4.12. Assume that the system (|6]) has a positive homogeneous steady state 
solution {u^,v^) and suppose that J and a > satisfy the following conditions: 

(1) There is an with 1 < i^, < N such that 
(a) Oj. > 0, and 
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(b) a directed path from Vi^ to Vn in G{J^) exists. 
(2) J IS Metzler. 

Then, if either the product xu* or is sufficiently large, (-u^j-u*) is linearly unstable. 

Proof. As in the proof of Theorem I4.8[ fix < 0, choose r large enough so that B = 
M{fi) + rl has positive diagonal entries, and consider the graph G{B^) with vertices 
Vo,Vi,...,Vn. 

We proceed by identifying a cycle in the graph. Since [-B*]o,i, = Oi. > 0, there is an edge 
from Vo to Vi^. Also, by our assumption, there is a path from V^, to Vat. Moreover, since 
[i?*]Aro = K > there is an edge from to Vq, thus completing a cycle. In particular, 
we have shown that Vq, Vi^, and are part of the same strongly connected component 
in G{B'). 

Therefore, finding a permutation matrix P as in ( !T7|) . we conclude that there is a 
triangular matrix T* = PB^P^ with a diagonal block Tj corresponding to the strongly 
connected component containing Vo,Vi^, and Vn. Due to the strong connectivity, Tj is 
irreducible. Now, we can prove that 

lim p(Tj) = oo and lim p(Tj) = oo, 

by repeating the arguments in the proof of f ll4p . because Tj is a nonnegative irreducible 
matrix with a structure similar to the matrix in (IT^ . 

Thus, for sufficiently large K or ctj^, the triangular matrix T has an eigenvalue larger 
than r, and so does B. Hence M{fi) = B — rl has a positive eigenvalue. □ 

5. Examples 

In this section, we present several examples of how to use the theory developed in the 
previous sections. We consider the system ([6]) with a g{v) that describes the kinetics of 
some CRN linking the chemicals wi, f2, . . . , f at. In all the examples, g is obtained by the 
laws of mass action kinetics. 

Example 5.1 ((The linear case) ). Suppose g{v) is linear, i.e., 

(18) g{v) = Av 

for some N x N matrix A. This limits the network to reactions of the form Vi vj ("vj 
produces f/'), or Vi — )■ 0, ("fj decays"), but these component reactions can be combined 
arbitrarily. 

Let us first observe a consequence of mass action kinetics. The reaction Vi — )■ Vj with 
a rate constant kji gives rise to an off-diagonal entry kji in the matrix A. Since the rate 
constants are nonnegative, this implies that A must be a Metzler matrix. 

If, in addition, we know that —A is a nonsingular M-matrix, then Proposition 13.41 can 
be applied to conclude the existence of nonnegative steady states. 

To investigate stability, we proceed by studying the paths in the CRN and applying 
our previous results, as exemplified next. 

Corollary 5.2. Assume that has a positive homogeneous steady state solution (m*, w*). 
Assume that the kinetics of a CRN on the chemicals Vi,V2, ...jVn is described by a linear 
function g{v) which is obtained by the law of mass action kinetics. If, for some chemi- 
cal Vi^, there is a path in the CRN from Vi^ to vjy, and at^ > 0, then {u^,v^) is linearly 
unstable whenever xu* or is sufficiently large. 
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Proof. We only need to verify tlie conditions of Theorem l4.12[ View tlie CRN as a directed 
grapli on vertices Vi, V2, . . . , Vat with an edge from Vi to Vj if and only if there is a reaction 
producing the chemical vj from the chemical Vi. This graph, except possibly for diagonal 
loops, is exactly G{A^), where A is as in ( lT8l) . The linearity of g also implies that the 
Jacobian J in (fTT]) coincides with A. Hence, the given path condition verifies Condition ([1]) 
of Theorem 14. 12[ We have already seen above that the assumptions of mass action kinetics 
imply that J = A is Metzler, thus verifying Condition (|2]) of Theorem 14.121 □ 

Example 5.3 ((Dimerization and decay) ). One of the simplest chemical reactions is 
dimerization, or the reaction 

(19) 2vi < — > V2. 

Additionally, suppose the chemicals Vi and V2 decay at rates 71 > and 72 > 0, re- 
spectively. Let the rate constant for the forward reaction 2vi — )■ V2 be ki > 0, and of 
the reverse reaction be ^2 > 0. Then, according to the law of mass action kinetics, the 
function g describing the CRN is 

-kivf + k2V2 - 7if 1 

kivj - k2V2 - 72^2 

Finally, assume that a species of density u produces the chemical vi at a rate of ai > 0, 
but does not produce V2 {0:2 = 0) and consider ([6]) in this setting. 

Recall from Section [3] that homogeneous steady states are solutions of ([8]), namely, 

(20a) aiU — kiv^ + ^2^2 ~ liVi = 0, 

(20b) kivf - {k2 + -f2)v2 = 0. 

Let us verify that this system admits positive steady states. Note that upon adding them, 
the equations of (!20l) are equivalent to: 

(21a) aiU = 7if 1 + 72f 2, 

(21b) V2 = TT^vl 

h + 72 

Eliminating V2 from (12 lap using (121b|) . we find that Vi satisfies 

fcl72 o 

22 r^^^i + 71^1 - am = 0, 

^2 + 72 

a quadratic equation which has a unique positive root whenever u > 0. Letting r{u) 
denote this positive root, we find the family of homogeneous stationary states for this 
example: 

ki 

u > (arbitrary), Vi = riu) > 0, f 2 = r'^iu) > 0. 

k2 + 72 

We next turn to determining stability for this example. Clearly g is not linear, so we 
cannot apply Corollary 15.21 However, we can apply Theorem 14.121 The Jacobian of g 
evaluated at a positive homogeneous steady state (m*,-!/*) is 

J ^ -2kiv^^i - 7i k2 

[ 2kiv^^i -k2 - 72_ 

where is the first component of v^. Then J is a Metzler matrix with a path in 
G{J^) from f 1 to V2 since 2kiv^:^i > 0. Hence Theorem 14.121 implies that any positive 
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homogeneous steady state of IQ, with the above g, can be destabihzed for large 

enough values of xu* or ai. 

Notice that the Metzler matrix J above has two negative eigenvalues (since the trace 
of J is negative, and the determinant is positive). Therefore, the CRN by itself (without 
the introduction of chemotaxis) is stable. However, once we introduce chemotaxis, the 
arguments above show that with a high enough chemotactic sensitivity (xu*), or a high 
enough chemical production by the species (ai), the stable system can be destabihzed. 



Example 5.4. Now we give an example where a hypothesis of Theorems 14.81 and 14.121 
fails. Nonetheless, we can analyze the system by applying Theorem 14.11 
Consider a reaction of the form 

Vi + V2 < > f 3 

where the rate constants of the foward and backward reactions are given by fci > and 
k2 > respectively. Again, assume that a species of concentration u produces only the 
chemical vi at a rate ai > and assume decay rates 7^ and production rates ai of chemical 
Vi satisfy 

7i > 0, 72 = 73 = and ai = a > 0, 02 = 03 = 0. 
Then the function g describing the reaction network is 



-kiVxV2 + fc2f3 
k\V\V2 - k2V3 



To verify that positive steady states exist, we look for u > and v > satisfying 

(23a) —kiViV2 + k2V3 — jiVi + au = 0, 

(23b) -kiViV2 + k2Vs = 0. 

Notice that the third equation of g{v) + au = is omitted because it is the same as the 
second up to the sign. By subtracting fl23al) from (123b|) we obtain the equivalent system 



(24a) 
(24b) 

From these equations it is clear that 



7ifi = au, 

klViV2 = k2V3. 



u > (arbitrary). 



a 



Vl 



-u, f 2 > (arbitrary). 



^^3 



aki 



-UV2- 



7i ' " 71^2 

is a positive solution of (123|1 . 

Having verified that positive steady states exist, we turn to determine their stability. 
The Jacobian of the function g evaluated at a positive homogeneous steady state {u^,v^) 
is 

-kiv^^2 - 7i -kiv^^i k2 

J = — fclf*,2 —klV^^i /C2 

/^lf*,2 -k2 

Since J is not Metzler (J12, J21 < 0), a hypothesis of Theorems I4.8l and l4.12l fails. 
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However, by virtue of Theorem 14 .11 the hnearized stabihty of steady states is determined 
by the eigenvalues of the matrix 



M(/x) 












71 





—kiv*,i 



-XU*IJ- 

-k2 + /il^s. 



-kiv^^2 + fJ'Di 
-kiv^^2 

klV^^2 

as /i varies over the spectrum of the Laplace operator. Fix any /i < and let K = —xu^,^. 

First consider the case K = 0, when there is no chemotactic term. Then, M(/x) is block 
triangular and its eigenvalues are fiD < and those of its next diagonal block, which can 
be expressed as 

-a — di —h c 
—a —h — d2 c 
a b —c — 



using the positive numbers a = fcif*,2, b = fcif c 
c?3 = —fiD^. Its characteristic equation is 



k2, di = 'ji — fiDi, d2 = — and 



+ b2y + bi\ + bo :-- 



+ {a + b + c + di + d2 + (i3)A^ 
+ {a{d2 + ds) + b{di + d^) + c{di + 6/2) + did2 + did^ + ^2^3) A 
+ {ad2d3 + bdid-^ + cdid2 + did2d^) = 

We apply the Routh-Hurwitz criterion [9], which says that all solutions of the character- 
istic equation have negative real part if and only if the following inequalities hold: 

62 > 0, 60 > and 6162 - 60 > 0. 

The first two are immediate. The third can be verified after some calculations: every 
term in the product 6162 is positive, and every term in 69 cancels some term (but not all 
terms) in 6162- Hence for any /i < all eigenvalues of M(/i) have negative real part. 

We next turn to the case where K > 0. Using cofactor expansion along the first row, 
we can compute det M{fi) to be 



C-Kdet 



ai -fcif*,2 + fJ-Di 

-kiv^^2 
kiv^o 



71 



—kiv^^i 
-kiv^^i + HD2 

klV^A 



where C = —fiDbo is the determinant of the matrix M{fi) in the case K 
C > is independent of K. Simplifying and identifying the K dependence, 

det M(/i) = C + Kaikiv^^2fi'D2 



0. Clearly 



where Kaikiv^,2fJ'D2 < 0. Hence there is a Kq > such that detM(yu) < for all 
K>Ko. 

To conclude, we claim that M(yu) must have a negative eigenvalue when K > Kq. If 
all eigenvalues of M(/x) are real, the claim is obvious. If not, complex eigenvalues occur 
in conjugate pairs, so in order that their product is negative, at least two of the four 
eigenvalues must be real. Furthermore, these two real eigenvalues must have opposite 
signs to satisfy det M(yu) < 0. 

Therefore, for all sufficiently large values of -fC, the steady state is linearly unstable. 
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